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Stability Results for Scattered Data 
Interpolation by Trigonometric Polynomials 

■ 
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A fast and reliable algorithm for the optimal interpolation of scattered 
. data on the torus T*^ by multivariate trigonometric polynomials is presented. 

I The algorithm is based on a variant of the conjugate gradient method in 

combination with the fast Fourier transforms for nonequispaced nodes. The 
main result is that under mild assumptions the total complexity for solving 
the interpolation problem at M arbitrary nodes is of order (D{MlogM). 
This result is obtained by the use of localised trigonometric kernels where 
the localisation is chosen in accordance to the spatial dimension d. Numerical 
^ , examples show the efficiency of the new algorithm. 
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i!^ ' 1 Introduction 

^ ' We discuss the approximation of scattered data by d-variate periodic functions / : T 

C, where T := [— |, ^) denotes the torus. In practical applications we are often con- 
fronted with the situation that experimental data or measured values of a function are 
^ ' only known at a finite sampling set X := {xj G T"^ : j = 0, . . . , M — 1}. Especially, 

■ nonuniform sampling sets appear in more and more applications in recent years. Given a 
notion of the distance of two points by dist {x, xq) := min^g^d || {x + j) — xq\\^, we mea- 
sure the "nonuniformity" of X by the mesh norm and the separation distance, defined 
by 

(5:=2max min distfa^i, a;), q := min dist (a;,-, a;/) , 
respectively. Obviously, the relation q < M^^^'^ < 5 is fulfilled. 
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For given samples {xj,yj) G T'' x C, j = 0, . . . , M — 1, a polynomial degree N G 2N, 

and the index set In := { — ^^■■■■>Y~^}'^ frequencies, we construct a d-variate 
trigonometric polynomial 

f{x) := J2 Ae'"'*'" 
fee/jv 

such that f{xj) yj, j = 0, . . . , M — 1. Turning this into matrix vector notation, we 
aim to solve the system of linear equations 

Af^y (1.1) 

for the unknown vector of Fourier coefficients / := {fk)keiN ^ C^'' . Throughout the 
paper, we denote the vector of the given sample values by y := {yj)j=o,...,M-i G and 
the nonequispaced Fourier matrix by 

\ / j=0,...,M-l;fce/iv 

In contrast to the widely used nonequispaced FFT for the fast matrix vector mul- 
tiplication with A, see [16] and its references, the efficient solution of (1.1) is still a 
challenging goal. Besides recently developed sparse reconstruction techniques, see e.g. 
[5. 13] and their references, a standard method to determine / is to solve the general 
linear least squares problem ||/||2 min subject to \\y — Af\\2 = min, see, e.g., [3, 
p. 15]. This can be done by means of the singular value decomposition which is not 
practical in the present situation for large problems due to its time and memory require- 
ments. Direct solvers for the univariate case d = 1 in [6, 18] obtain a solution in 0{NM) 
floating point operations. 

For N"^ < M, the linear system (1.1) is over-determined, so that in general the given 
data y will be only approximated up to a residual r :=y — Af. In order to compensate 
for clusters in the sampling set X, it is useful to incorporate weights Wj > and to 
consider the weighted approximation problem 

M-l 

\\y - Affw = E ""j^yj - ^ (1-2) 

j=0 

where W := diag(w;j)j=o,...,M-i- In [10] it has been proven that this problem has a 
unique solution if iV < {]^dS)~^. Its solution is computed iteratively by means of 
the conjugate gradient method in [7, 2, 9], where the multilevel Tocplitz structure of 
A'^VFA is used for fast matrix vector multiplications. Slightly more stable with respect 
to rounding errors is the CGNR method, cf. [3, pp. 288], which iterates the original 
residual ri = y — Afi instead of the residual A^Wri of the normal equations. Note 
furthermore, that it has been suggested in [17] to incorporate some "knowledge on the 
decay of the Fourier coefficients" W := diag{wk)keiN^ ^fe > 0. Their approach is based 
on the weighted least squares problem 

\A''w{y-Af^ 
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In contrast, we focus on the under-determined and consistent linear system Af = y, 
i.e., we expect to interpolate the given data yj G C, j = 0, ... ,M — 1, exactly. We 
show that the nonequispaced Fourier matrix A has full rank M for every polynomial 
degree N > 2dq~^. In particular, we incorporate damping factors Wk > 0, fe G In, and 
consider the optimal interpolation problem 

Wffw-^ = V ^ ^ min subject to Af = y, (1.3) 

where W := diag(t(;fe)feg/^. We prove that for a large class of "smooth" damping factors 
Wk problem (1.3) is well-conditioned, where the "smoothness" has to be chosen with 
respect to the spatial dimension d. We propose to solve problem (1.3) by a version of 
the conjugate gradient method in combination with the nonequispaced FFT [19, 16, 12] 
to efficiently perform each iteration step. 

The outline of this paper is as follows: In Section 2 we set up the basic notation 
and relate the optimal interpolation problem (1.3) to a particular trigonometric kernel. 
Furthermore, wc propose Algorithm 1 for computing the solution of the interpolation 
problem efficiently. For the sake of analysing the convergence of this algorithm, we then 
present our theory on localised trigonometric kernels in Section 3. Our first result in 
Theorem 3.3 is a version of the typical smoothness-decay principle in Fourier analysis 
and relates the "smoothness" of the weights in (1.3) to the localisation of the correspond- 
ing trigonometric kernel. We use this decay in Section 4 to prove that well separated 
sampling nodes yield a stable interpolation problem (1.3). The eigenvalue estimates 
are given for the univariate setting in Theorem 4.1 and for the multivariate setting in 
Theorem 4.6. Subsequently, Corollary 4.7 applies the general result to a particular class 
of damping factors and concludes with conditions sufficient for the full rank of A. As 
the equidistant case in Theorem 4.10 and Corollary 4.11 reveals, the assumption on the 
separation distance is of optimal order. We provide numerical examples in Section 5 and 
draw our conclusion in Section 6. The software and all numerical examples are available 
from our NFFT- homepage [12]. 



2 Optimal interpolation and its iterative solution 

After setting up our notation in Definition 2.1, we prove in Lemma 2.2 that the optimal 
interpolation problem (1.3) can be stated as normal equations and the matrix in these 
equations obeys special structure. Furthermore, we propose Algorithm 1 for the iterative 
solution of the interpolation problem and state a basic convergence result for this scheme. 

Definition 2.1. Let d G N, e 2N, and In = {-f , . . . , f - l}'' be given. We define 
for positive weights Wk > 0, fee In, with normalisation ^keiN ^fe ~ -'^ ^^'^ x ef^ 
the trigonometric kernel 

Kn{x) := ^fce2-*=-. 

fcG/jv 
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The particular class of tensor product kernels is given by 



d-i 

t=0 

where Kn denotes a univariate kernel and x = (xq, • • • , Xd-i)~^ ■ 

Furthermore, given a sampling set X cT'^, we define the kernel matrix 

Kn := {Kn{xj - a^;)),. ,=o,...,M-i ^ C^"^- (2-1) 

We denote by A = A (Kn) and A = A (Kn) the largest and smallest eigenvalue of 
the kernel matrix Kjsi, respectively. Their ratio is denoted by the condition number 
cond(Kjv) = j. □ 

Note, that from the definition immediately follows Ki\f(0) = max^^fd \K]\f{x)\ = 1 
and (K]sf)jj = 1, j = 0, . . . , M — 1. The following theorem collects some basic facts. 

Lemma 2.2. Let the number of nodes M € N, the sampling set X C T'^, the polynomial 
degree N G 2N, and the damping factors Wk > 0, k & In, be given. The optimal 
interpolation problem (1.3) is equivalent to the damped normal equations of second 
kind 

KNf = y, f = WA^'f, (2.2) 
where the kernel matrix Kn & ^m^m Q^gyg ^|jg factorisation 

Kn = AWA^, (2.3) 

hence is positive semidefinite. 

Proof. The second assertion follows from {AWA^)j^i = J2kelN e^'^'*'"'^ WfcC^^'^''"'^' and 
Definition 2.1. Furthermore, a solution / of Af = y has minimal weighted norm if 
and only if it is perpendicular with respect to the weights to the null-space of A, i.e., 

- —1/2 - - 1/2 

W / _L J\f{AW ). We conclude (2.2) by the fact that the orthogonal complement 
of the null-space of a matrix is just the range of its adjoint. ■ 



Denoted in Algorithm 1 by CGNE, cf. [3, pp. 288], we solve the Normal equations 
(2.2) by the Conjugate Gradient method, minimising in each iteration the Error. 

The proposed method finds approximations from a Krylov space closely related to the 
one of the CGNR method for (1.2), but with minimal error instead of minimal residual. 
Note that we exploit the factorisation in (2.3) to iterate the original vector / instead of 
the vector /, cf. equation (2.2). Hence, we use fast matrix vector multiplications for A 
and A^ by means of the fast Fourier transforms at nonequispaccd nodes (NFFT) having 
an arithmetical complexity of 0{N'^ log(N'^) + M\ loge|^) in each iteration, where e is 
the prescribed accuracy. Details concerning NFFT algorithms can be found for example 
in [19, 16] and a corresponding software package in [12]. Applying the standard estimate 
for the convergence of the conjugate gradient method we obtain the following lemma. 
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Algorithm 1 CGNE 

Input: dimension d G N, number of samples M G N, polynomial degree iV G 2N; 
sampling set C T'^, samples y G C*^, and initial vector /q G C'^'' 



ro = y - A/o 

for Z = 0, ... do 

fl+i = fl + aiWpi 
ri+i = ri- aiAWpi 

Pi+i = PiPi + A^n+i 

end for 

Output: the Z-th iterate fi 



Lemma 2.3. Let the kernel matrix K n in (2.1) be regular and let ei := fi — WA^Kj^^y 
denote the error of the l-th iterate within Algorithm 1. Then the a-priori error bound 

holds true. 

Proof. We note that ||e;||^-i = ||/; — Kj^^y\\Kf^, where fi denotes the Z-th iterate 
of the conjugate gradient method applied to equation Kj^f = y, cf. Lemma 2.2, and 
apply the standard estimate for the conjugate gradient method, see also [3, pp. 288]. ■ 

This result includes the special case of M = A^^^ equispaced nodes and no damping, 
i.e. = 1, fc G In, where the first iterate of our algorithm is already the solution to 
equation (2.2). We present estimates for the extremal eigenvalues A, A dependent only 
on the quantities N,q, and the damping factors Wk, k £ Analogous results for 

the stability of the interpolation by radial and zonal functions are obtained in [15, 20]. 
Section 3 prepares our estimates by constructing localised kernels. 

Remark 2.4. Before that, we would like to comment on the following: 

The weighted norm in (1.3) is induced by the inner product g^W f. In particular, 
the definition {f,g)~-i := g^W f makes the space of trigonometric polynomials 



5 



Tat := span I e^'^'*'' : k G /tv} to a reproducing kernel Hilbert space. Its reproducing 
kernel is given by Kj^, i.e., the point evaluations obey f{x) = {f,K]y[- — x))^-i. 

Moreover, the solution f{x) = J2keiN /fc^^'^'*"" normal equations (2.2) has 

comparable norm to the given samples, i.e.. 

Ml < if, g)^-i<\-' Ml 

This norm equivalence is due to fact, that the field of values of the matrix KJ^^ is 
bounded by its extremal eigenvalues and y^Kj^ y = f K^f = f W f. □ 



3 Localised kernels 

Starting from a class of admissible weight functions in Definition 3.1, we construct lo- 
calised trigonometric kernels in Theorem 3.3, where Lemma 3.2 serves as an interme- 
diate step. Following the smoothness-decay principle in Fourier analysis, we relate the 
smoothness of the weight function to the decay of the kernel Kn built upon the sampled 
weights. A related approach is taken in [14, Thm. 2.2] for the detection of singularities. 
The particular class of B-Spline kernels, cf. Definition 3.4, is considered in Corollary 3.5. 
While we present our results on the connection between smooth weight functions and 
localised kernels for the univariate case, we give its generalisation to the class of tensor 
product kernels in Corollary 3.7. 

Definition 3.1. For /3 G N, /? > 2, a continuous function : M ^ M is an admissible 
weight function of order (3 if it is nonnegative, possesses a (/3 — l)-fold derivative g^^~^^ 
of bounded variation, i.e.. 



V 



/n—1 
d^^^-i) (^)|=sup5^ 



■j+i. 



< oo, 



where the supremum is taken over all strictly increasing real sequences {zjjjgNj,, and 
satisfies the additional properties supp^f = [— ^, g^'^\±^) = for 7 = 0, . . . , /3 — 1, 



< ^, and the normalisation 



\g\\Li = 1. We denote by BV^ ^ the set of 



g{z) > 0, \z 

admissible weight functions of order /3. 

Furthermore, we define for notational convenience the zeta function ^(/?) := ^ 
/3 > 1, and for g G BVq ^ the norm of the samples 



00 

r=l 



li,Ar 



□ 



The following lemma prepares Theorem 3.3. 
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Lemma 3.2. For /3 G N, (3 > 2, let a function g G SV^j'^ ^ be given. Then for N G 2N, 
N > 2/3, and x G [—5,5] \ {0} the following estimates hold true 



N_ 

2 



2mkx 



< 



(2/3-l)C(/3)|5(^-^)| 



V 



li,iv > iv(l-2C(/3)(4vr/3)-^|5(^-i)|J 



Proof. First, wc define for 2; G [— the function hx{z) := g (z) e^'^^^^^ . Thus, 
the Poisson summation formula yields 



2 



(z) e-^^^'^^'^dz 



and by applying integration by parts and the fact that g^'^'> (±5) = for 7 = 0, . . . , /3 — 2 
further 

f ^ 

^ /i^ = AT ^ (2vri.V (r - x))-^'^-^) / g^^-'^ {z) e^-Ar.C.-.)^^ 

^ 2 



(27riAr)'^ 



1 

2 

^ (r - x)-'^ / 5^'^"^) (^) e^-^^^^^--^)^ d. 
r-ez _i ^ ^ ^ 

2 



< 



r-ez 

1 + \x\^ ^ |r — xp'^ 

r€Z\{0} 



a sup 

(27r)^7V/3-i|x|/3 rosz 



1 

2 

_ 1 
2 



Using 1 + |x|^E^g2\|o| V - < (2^ - 1)2^"^C(/?) for |x| < i and 



dz 



dz 



< 



yields the assertion. 

By the Poisson summation formula, we note furthermore that 



^Il5lli,jv>l 



-2iTiNrz 



dz 



reZ\{0}_i 



and proceed analogously in order to prove the second assertion where we use > 2/3 to 
obtain an estimate independent of N. ■ 
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Theorem 3.3. For /3 G N, (3 > 2, let a function g e BVq ^ be given. Furthermore, let 

N G 2N, N > 2(5, and the damping factors 



he given. Then the kernel Ki^, cf. Definition 2.1, fulfils 
\Kn{x)\ < 



N N 



(2/3-l)C(/3) 




1) 


V 


1 


2/3-1 (27r>^-C(/3)/?-'3 


y(/3-l) 





forxe \{0}. 
Proof. Note first, that 



Kn{x) 



1 + e- 



N_ 

-2mx 2 



2N 



E 



N 



Jlmkx 



Thus, we obtain i^Ar(O) = 1 and by applying Lemma 3.2 also the decay property. 

We apply Theorem 3.3 in the following to the particular class of B-Spline kernels. 
Definition 3.4. Let /? G N be given. The normalised B-Spline is defined by 



gi3 (z) := PNp Pz + 



where Np denotes the cardinal B-Spline of order /3. The cardinal B-Splines arc given by 
Ni{z) = 1 for z G (0, 1), Ni{z) = elsewhere, and N(s+i{z) = Np{t)(1t, see e.g. [4]. 
Furthermore, we define for /3 G N and AT G 2N the B-Spline kernel by 



Bp,N {x) ■= 



1 + e 



JV 

-2mx 2 



1,7V 



□ 



Corollary 3.5. Let /3 G N, [3 > 2, and iV G 2N, > 2/3, be given. Then the B-Spline 
kernel Bp^^^x), cf. Definition 3.4, fulfils 

(2.«-l)C(/?)/3^ 



\Bp,N {x)\ < 



\Nx\- 



2/3-i7r/3-C(/?) 
for X G [-5, |] \ {0} and Bg^iO) = 1. 

Proof. Note that gp G BV^'K Using N'p{z) = Np_i{z) - Np^i{z - 1), we conclude 

T=0 ^ 



V 



V 



(5^ 



and apply Theorem 3.3. 
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Note, that in contrast to [15] the order (3 of the B-Spline and the degree N G 2N of 
the kernel -B^,Ar are independent of each other. The special case /3 = 1, i.e. the "top- 
hat function" gi{z) = 1 for \z\ < ^ and gi{z) = elsewhere, leads to the well known 
Dirichlet kernel Bi^]\[{x) = jj ^^.g/^^ ^^rnkx^ Analogously, [3 = 2, i.e. the "hat function" 
52 (^) = 2 — 4|z| for < ^ and g2{z) = elsewhere, leads to the Fejer kernel. The 
increasing localisation of the B-Spline kernels is illustrated in Figure 3.1. 




-0.5 0.5 -0.5 0.5 -0.5 0.5 



Figure 3.1: From left to right: Real part of the Dirichlet kernel Bi^20, the Fejer kernel 
-62,205 f'nd the B-Spline kernel 54^20- 



Remark 3.6. If we assume in Corollary 3.5 furthermore, that A'^ = Pa, a G N, then the 
stronger estimate \Bfj^N{x)\ < 2C(/?)(1 — 2~^)(^)^|A''x|~^ holds true. This improvement 
is due to 11^/3 II i,Ar = N in Lemma 3.2 and follows from the partition of unity of the 
cardinal B-Spline Np and the refinement equation Nfj{z) = Y1t€1^^^''^'' -^pi^^ ~ ''") 
some finitely supported coefficients a^'^''^^ > 0, see e.g. [4, pp. 8]. 

In particular, the Fejer kernel fulfils \B2,n{x)\ < |Ax|~^, N G 2N, which also follows 
from the estimate |sin(7rx)| > 2|x| for |x| < 1/2 and the representation 

2(14- e-2^'^) / sin (f vrx) \ ^ 



B2,N (x) = 



N'^ \ sin (ttx) 



Along the same line follows the localisation property Bi^n{x) < \Nx\ ^ and Bi^jv(O) = 
1 for the Dirichlet kernel. □ 

We complete this section by an extension of our result to the multivariate case d > 1. 
Indeed, tensor products of the kernels constructed in Theorem 3.3 yield also localised 
multivariate kernels as shown in the following corollary. 

Corollary 3.7. Let the univariate kernel Kn, cf. Definition 2.1, fulfil for some /3 G N, 
some constant > 1, and x G [— ^, 5] \ {0} the decay condition \Kn{x)\ < C^|A/"x|~^, 

then its tensor product kernel Kn{x) = nf=o ^N{xt) fulfils for x G [— ^, \ {0} the 
estimate 

' ^ ^' " A^^lla^ll^ 

Proof. The assertion follows simply from the estimate \Kn{x)\ < |.RrAr(||a;||oo)|- ■ 
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4 Stability of the interpolation problem 



The noneqTiispaccd Fourier matrix A has fTill rank M for d = 1 and N > M. Unfor- 
tunately and due to the famous result of Mairhuber-Curtis, see e.g. [20, Thm. 2.3], we 
cannot expect full rank of the matrix A for N > Md in the multivariate case d > 1. 

However, we prove stability results for the trigonometric interpolation problem at q- 
separated nodes in the univariate case, cf. Theorem 4.1 and the multivariate case, cf. 
Theorem 4.6. Basically, a localised kernel Kj^ yields a diagonal dominated kernel matrix 
Kn and thus full rank of the nonequispaced Fourier matrix A for g-separated nodes. 
Furthermore, we prove stability results for a slightly generalised interpolation problem 
at equispaced nodes and subsets of equispaced nodes in Theorem 4.10. These results are 
applied to the B-Spline kernels from Section 3 in Corollary 4.7 and Corollary 4.11. 



The univariate setting 

The following Theorem 4.1 gives estimates for the extremal eigenvalues of the matrix 
Kn' under reasonable assumptions on the kernel K]\f. 

Theorem 4.1. Let N G 2N be given and let the kernel Kn, cf. Definition 2.1, fulfil for 
some (3 > 1 and x e [—^i^] \ {0} localisation property 

\Kn{x)\< 



Nl^\x\P ■ 

Furthermore, let a sampling set X contain arbitrary nodes with separation distance 
q> 0. Then, the extremal eigenvalues of the matrix Kn are bounded by 

1_?<(|)££<A<1<A<1 + ?<(|1^£. 

Nf^qP - - - - NPqf^ 

Proof. As usual, let M denote the number of nodes in X. Due to (0) = 1, cf. 
Definition 2.1, we obtain trace(-fCAr) := X^^q^ Kn{0) = M. Since the trace is invariant 
under similarity transforms, all eigenvalues sum up to M and thus, the inequality A < 
1 < A is fulfilled. 

Now, let A* be an arbitrary eigenvalue of K]\f, then for some index j G {0, . . . , M — 1} 
Gershgorin's circle theorem yields 

M-l 

|A*-1|< \KNixj-xi)\. 

By using that the separation distance of the sampling set is q and by the localisation of 
the kernel Kn, we obtain 

W M ^ '^f^ 1 ^ '^Cfs ,-0 ^ 2C(/3)C/3 



^ \x - -x,f ~ Nf^qf^ ^ N^qli 

l=0-l^j \-^3 ^ 1=1 ^ 
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Indeed, the kernels constructed in Theorem 3.3 yield well conditioned matrices Kj^. 
We also note that the decay of the Dirichlet kernel only allows for a weaker result. 

Remark 4.2. Let a sampling set A' € T with separation distance q > he given. 
Then the application of Theorem 4.1, where the last step of its proof is replaced by the 
estimate Yll^/"^^ < 1 + ln^, yields: The matrix {Bi^p^{xj — x/)j,;=o,...,M-i = -k^-^ 
is nonsingular for A/" > (1 + \\og{2q)\) q"^ . The logarithmic term in this condition is 
clearly suboptimal. □ 

As an immediate consequence of Theorem 4.1 we state a stability result for an equis- 
paced grid disturbed by jitter. 

Corollary 4.3. Let the assumptions of Theorem 4.1 hold true. Furthermore, let the 

sampling nodes be of the form a^j = — | + "^"^j J = 0, . . . , M — 1, where < £j < £ < 1. 
Then the eigenvalues of the matrix are bounded by 



Proof. Since the separation distance is bounded hy q> M ^(1 — e) the result follows 



The multivariate setting 

First, we borrow a packing argument on the sphere from [15] and refine it in Lemma 
4.5 for the present setting, i.e., we show how many g-separated nodes can be placed in 
a certain distance to a reference node, see also Figure 4.1. 



1- 



2(:{P)CpM^ 
NP (1 - ef 



<A<1<A<1+ 



NP (1 - ef 



by Theorem 4.1. 




for m = 0, . . . , Lg~V2j - 1 and 



^9,U-V2J ■= {a; e T'^ : [g-V2j « < dist (a;,0) < I/2} . 



Its restriction to the sampling set X will be denoted by Rx 



,q,m 



■.= R 



■q,m 



□ 



Lemma 4.5. Let d G N and an q-separated sampling set X with < q < ^ be given. 
Then, each of the sets Rx,q,m ^Sts hounded cardinahty 



\Rx,q,m\ < (2'^ - 1) m'^'\ m = 1, . . . , [q-^/2\ . 



11 




r 



L 




i 

■J 



Fi gure 4.1: Partitioning of the torus into the rings Rq^rm = 0, . . . ^/2J (left). 

Further subdivision into shifted and rotated versions of the cube [0, g)*^, where 
arrows indicate the "ownership" of the faces to a particular cube (right) . 

Proof. We use a packing argument for the partition {-Rg,m) m = 0, . . . , [q'~^/2J} of 
the torus T*^. Each ring Rq^rn is subdivided into shifted and rotated versions of the cube 
[0, g)'^, of. Figure 4.1 (right). This is done such that each point in Rq^m is contained in 
at least one of these boxes and the boxes share no interior points with each other. Every 
box contains at most one node of the sampling set and hence, the estimate 

\R?^,,,m\<^, I da;<2'^((m + l)^-m'^)=2''^Qm'^-*<2'^m'^-iX:(^)- 

is valid. ■ 

Using localised kernels in conjunction with a separated sampling set and Lemma 4.5, 
we state the following theorem on the stability of the interpolation problem. 

Theorem 4.6. Let (i G N, iV G 2N be given and let the kernel Km, cf. DeGnition 2.1, 
fulSl for some (5 > d and x & {0} the localisation property 

\Kn{x)\< 



Furthermore, let a sampling set X contain arbitrary nodes with separation distance 
< < 5. Then, the extremal eigenvalues of the matrix Kn are bounded by 

2d (2^^ - 1) a/? - + 1) Cg 2<^ (2^^ - 1) C (/? - + 1) Cr 
1 5^ '-^^^ ^ < A < 1< A < 1 + ^ '-^^^ 

Proof. Let A* be an arbitrary eigenvalue of K^- Without loss of generality, let the 
diagonal element of the matrix used in Gershgorin's circle theorem correspond to 



12 



Xo = 0. Then we conclude by -fCjv(O) = 1, cf. Definition 2.1, that 



M-l 



|A,-1| < ^ |i^iv(0-xO|. 



Using the partition from Definition 4.4, Lemma 4.5, and the localisation of the kernel 
Kn, we get 

m=l XieRx,q,m. 

2^ (2^ - 1) Cfs ^'^'^ ,1 „ „ . 

m=l 

2'i (2'^-l)C(/3-d+l)C^ 

Particularly, this result includes Theorem 4.1 if we set d = 1. 

Corollary 4.7. Let the dimension d G N, an arbitrary sampling set A" G T*^ with 
separation distance < q < ^, and a polynomial degree N € 2N, N > 2dq~^, be given. 
Then the nonequispaced Fourier matrix A has full rank. Moreover, the eigenvalues of 
the kernel matrix Kn = AWA^ obtained from the B-Spline kernel of order (3 = d + 1 
are bounded by 



0<l-l5^) <A<l<A<l+(-j . 

Proof. Note first that > 2/3. We apply Theorem 4.6 where we use the estimates for 
Cp given in Corollary 3.5 and simplify the involved constant. Hence, the full rank of A 
follows. ■ 



Thus, we have shown that the optimal trigonometric interpolation problem at q- 
separated nodes in d dimensions obeys a uniformly bounded condition number for a 
polynomial degree N > 2dq~^ and appropriate damping factors. The dependence on 
q~^ is optimal as the subsequent analysis of the equispaced case shows. However, the 
constant 2d is not optimal for high spatial dimensions. As pointed out in [2] for the 
related approximation problem (1.2), it is an open problem to improve on this. 

In summary. Lemma 2.3 and Corollary 4.7 assure in our situation a prescribed reduc- 
tion of the error ||e/||^-i in a constant number of iterations. Hence, if we assume an 

additional uniformity condition q = cM~d for the sampling set X, the total arithmeti- 
cal complexity of Algorithm 1 for solving (1.3) up to a prescribed error is bounded by 
O(MlogM). 
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Remark 4.8. Other frequently applied kernels also possess specific localisation proper- 
ties and thus, yield stable interpolation at g-separated nodes by means of Theorem 4.6 
as follows: 

We define for /3 G 2N, a G N, and iV = /3(c7 - 1) + 2 the Jackson kernel by 

2nix 



, , , l + e-^'^'^ {smiaTTx)\ 
2af^ \ sm [TTX) ) 



Being a normalised power of the Fejer kernel Bi^ia^ the coefficients Wfc of the Jackson 
kernel can be obtained by an iterated discrete convolution of the coefficients of the 
Fejer kernel, see [1] for details. In contrast, the B-Spline kernel relies on a continuous 
convolution. The Jackson kernel is localised as | J^^Ar(a;)| < {^^^\Nx\~^ for x G ) ^] \ 
{0} and fulfils J/3,7v (0) = 1. Hence, the tensor product Jackson kernel of appropriate 
order yields for AT G 2N, AT > 2.\dq-^ , the nonsingular kernel matrix 



( Jprd+l-l Y — jg/)) 

\ ^\ 2 l'-*^ J j,l=0,...,M-^ 



Secondly, it is well known that the weight 1 + {2iikf°' is associated to the squared 
Sobolev norm ||/||2 + H/'^'^^lll- For /3 G N and q;,7 > 0, a regularised and slightly 
generalised weight is given by 

^ (i ~ ^ ) 



7 + kl 

for \z\ < ^ and ga,(3,'y{z) = elsewhere, where the constant c^^^^-y is chosen such that 
lba,/3,7l|Li = 1- Here, the denominator generalises the weight 1 + (27rfc)^" and the 
nominator ensures 

5a,/3,7 e BV^ \ We define for G 2N the Sobolev kernel by 

iV 

/ h \ 

2-iTikx 



^a,l3,^,N [X) := — II 2^ ga,/3,^ — 

^\\9a,l3,l\\l,N , N V^/ 



e 



The kernel is localised as \Sa,i3,^,N {x)\ < Ca,/3,^\Nx\-^ for X G [-^ k] \ W and some 
constant c^^^^^ > and fulfils Sa,/3,'y,N (0) = 1. □ 

Results for equispaced nodes 

In the case of equispaced nodes we employ the fact that the matrix Kn is circulant. We 
present a slightly generalised result in the following Theorem 4.10. 

Definition 4.9. We define for d, n G N and weights it)fe G M, k e with Ylke'Z^ l^^l 
oo, the kernel 

K {x) := ^ Wfce^'^*'" 
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and by evaluating at the equispaced sampling nodes j = {jo, ■ ■ ■ ,jd-i) ^ In, the matrix 

□ 



Theorem 4.10. The matrix K in Definition 4.9 possesses the following properties. Its 
eigenvalues are given by 

for s & In- For tensor product weights Wk = Ylt=o ''^h' X^fcez < ^^^^ simplifies 
to 

d-l 

As (K) =n'^YlYl '^^t+nrt 

for s e In- Moreover, the extremal eigenvalues of 



Kr--={K 



^ ' ' j,ier 



are bounded by the extremal eigenvalues of K for any V d In- 

Proof. The matrix K is multilevel circulant and thus diagonalised by the Fourier 
matrix F„ = (e^'''*'^'/'*)j,fee/„. We calculate 



3,iei„ \ ^ J 



,27ritZ/n 



^ U)fe ^ g-27rii(s-fe)/n ^ ^27ri«(t-fe)/n 



ioT s,t e In and use 



p27rij(s-fe)/n ^ J ^'^ if ^ ^ ^'^^ 



e 



n 

otherwise. 



See also [15, Cor. 3.10, Thm. 3.11] for the univariate case. 

The second assertion is due to the Kronecker product structure of the matrix K in the 
case of tensor product kernels. The last assertion follows from the fact that removing 
a node is nothing else than removing its corresponding row and column in K and from 
the interlacing property for eigenvalues, see [11, pp. 185]. ■ 

Now, let d,n E N, N ^ 2N, and the equispaced sampling set X = i/„ C T'^ be given. 
A simple consequence of Theorem 4.10 is the fact that the kernel matrix Kjy = jjAA^ 
is singular whenever N < = n. Hence, the condition > 2dq~^ for the full rank 
of A is optimal with respect to q. Nevertheless, we also apply Theorem 4.10 to obtain 
positive results in the equispaced setting for the Dirichlet and the Fejer kernel. 
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Corollary 4.11. Let the dimension d £ N, a sampling set X with M > 2 equispaced 
nodes, and a polynomial degree N £2N with N > q-^ = Md be given. Then, for the 
Dirichlet kernel Bi^n and its tensor product version for d > 1, the extremal eigenvalues 
of the kernel matrix obey 

d / r AT 1 \ d 



([MX 

\ Nq J 



A < 1< A 



Nq J 



Furthermore, the Fejer kernel B2^n yields 

d 



1 



1 



<A<1<A< 1+ 



where equality holds for the outmost inequalities if A'^ = (2a + 1) q ^, cr G N. In partic- 
ular, the kernel matrix Ki\f is nonsingular. 

Proof. Throughout this proof, let n = q~^ and the damping factors be extended 
by «)fc = for k ^ In- We apply for cZ = 1 the first statement of Theorem 4.10 to 
the weights Wk = jf, k e In, of the Dirichlet kernel for N e 2N and to the weights 

Wk = -^(1 — ^^^"^^ ), k G In, of the Fejer kernel, if A" = {2a + 1) q~^, a £N, respectively. 

The assertion is little more delicate for the univariate Fejer kernel and A^ 7^ {2a + 1) q~^. 
We use the representation 



\B2,n{x)\ = J^Y1 ^' 



r=0 k= 



Now, let be an arbitrary eigenvalue of the kernel matrix Kn, then Gershgorin's circle 
theorem yields 



n-l 



USE 



B- 



2,N 



An 
iV2 



E 2 



r=0 



+ 1 



Since for Q := [^^^:^\ and R := ^ — ^ — the identity 



Q-l (.s + l)ri-l 



nQ+R 



r=0 

holds, we proceed 
-1 



ES-E E ^+EQ- 



s=0 r=sn 



{N-nf-{2 {R+l) 
8n 



n 



4n 

iV2 



r=0 





r 








-fi- 



+ 1 



< 



4n 

iV2 

Ar2 

iV2- 



{n-Nf -{2 {R+l)-n) 
8n 



2(i?+ 1) -n 
TV 



A^ 



The case d > 1 is due to the second statement in Theorem. 
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5 Numerical results 

In this section, we exemplify our findings on tlic stability of the optimal interpolation 
problem (1.3) and its iterative solution by Algorithm 1. 

The estimates for the condition number of the kernel matrix Kjq for equispaced nodes, 
of. Corollary 4.11, are shown in Figure 5.1 (left). For Nq G N and the Dirichlet kernel 
Si^TV the matrix Kj^is just the identity. However, using the better localised Fejer kernel 
-B2,iv improves the condition number already for N > \/?,q^^ when Nq N. 

We present the effect on the stability of the interpolation problem when the equispaced 
nodes are perturbed by jitter error, of. Corollary 4.3, in Figure 5.1 (right). We choose 
different sampling sets of size M = 1, . . . , 100 with equispaced nodes disturbed by 10% 
jitter error and evaluate the maximum condition number over 100 reruns for the Dirichlet 
kernel -Bi,6M and the Fejer kernel -B2,6M, respectively. The Fejer kernel produces a lower 
condition number which is also validated by the shown upper bound. These results 
confirm the theoretical results of Corollary 4.11 and Corollary 4.3. 




100 200 600 20 40 60 80 100 



Figure 5.1: Condition number of the kernel matrix Kjsi. Left: Condition number with 
respect to polynomial degree N = 100, ...,600, no weights, i.e., Dirichlet 
kernel (dash-dot); weight function g2, i.e., Fejer kernel (solid), and the es- 
timate of Corollary 4.11 (dashed); here, the number of equispaced nodes is 
M = 100. Right: Condition number with respect to the number of nodes 
M = 1, . . . , 100, the nodes are equispaced perturbed by £rei. = 0.1 jitter er- 
ror, the polynomial degree is = 6M; no weights, i.e., Dirichlet kernel ( + ); 
weight function g2, i.e., Fejer kernel ( x ), and its estimate by Corollary 4.3 
(dashed) . 

Furthermore, we apply Algorithm 1 using the NFFT software package [12] to recon- 
struct a univariate signal from randomly scattered data in Figure 5.2 and show in Figure 
5.3 the reconstruction of a bivariate signal from a glacier data set [8]. The main tool 
in our iterative algorithms is the NFFT, i.e., the fast matrix times vector multiplication 
with A and A^, respectively. Details concerning NFFT algorithms can be found for 
example in [16] and a software package can be found in [12]. 
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The reconstruction of the randomly sampled univariate signal shows the decay rates 
of our iterative scheme. The sampling set consists of M = 100 nodes separated by 
q = 4x 10"^ and we reconstruct with a polynomial degree N = 1000 and the Dirichlet-, 
Fejer-, B-Spline-, and Sobolev kernel. All schemes converge within 15 iteration where 
this is justified only for the Fejer- and the B-Spline kernel. 



10" 



10" 



10" 



X K X X X K 



10 



15 



10" 



10" 



10 



10 



15 



10 



10"' 



10" 



xx^xxxxxxxxxxx 



10 15 




Figure 5.2: Native error \\fi — WA^K]^^y\\^-i for the univariate interpolation problem 
with respect to the current iteration /. The number of samples is M = 100, 
the number of computed Fourier coefficients is N = 1000, and the separation 
distance of the nodes is g = 4 x 10~^. Top left: no weights, i.e., Dirichlet 
kernel; Top right: weight function g2, i.e., Fejer kernel, predicted decay rate 
(dashed); Bottom left: weight function g^, i.e., B-Spline kernel, predicted 
decay rate (dashed); Bottom right: weight function 51,2,10-25 i-e., Sobolev 
kernel. 

The last example shows a typical test case known in radial basis function methods. 
We reconstruct from a data set of M = 8345 samples on level curves of a glacier a total 
number of 2^ x 2^ ~ 8M Fourier coefficients. Note however, that the sampling set is 
highly nonuniform in the sense that the separation distance is very small compared to 
the mesh norm. The assumptions of Theorem 4.6 are not fulfilled. Nevertheless, the 
proposed method yields a very good approximation to the given data after 40 iterations, 
which is also supported by the cross validation test in Table 5.1. Here, we exclude 
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M randomly chosen samples X d X , y ^ ^ C from the reconstruction process and 
compare our approximations at these left out nodes. The comparison is done by means 
of the relative data residual and the relative validation residual after 40 iterations 

\\yh ' ' \\yh 

As can be readily seen, the CGNE scheme achieves both a small data residual r and a 
small validation residual f. The proposed CGNE method combines the good data fit 
of the CGNR scheme {N = 256) with the smooth approximation of the CGNR scheme 
(AT = 64). 




Figure 5.3: Reconstruction of the glacier data set vol87.dat from [8], M = 8345 nodes, 
N = 256, 40 iterations, tensor product damping factors Wk to the weight 
function (71 3^0-3; sec glacier in [12]. Left: surface plot. Right: contour 
plot and sampling set (•). 





CGNE 


CGNR (A^ = 256) 


CGNR {N = 64) 


M 


r 


f 


r 


f 


r 


f 


200 


6.9e - 04 


1.7e - 02 


5.0e - 04 


lAe - 01 


8.3e - 03 


1.7e - 02 


400 


4.7e - 04 


2.3e - 02 


5.0e - 04 


2.0e - 01 


8.3e - 03 


2.3e - 02 


600 


5.7e - 04 


2.9e - 02 


5.1e-04 


2.5e-01 


8.1e - 03 


2.9e - 02 


800 


4.7e - 04 


3.4e - 02 


5.0e - 04 


2.8e - 01 


8.0e - 03 


3.4e - 02 


1000 


4.6e - 04 


3.8e - 02 


4.7e - 04 


3.2e -01 


8.0e - 03 


3.8e - 02 



Table 5.1: Cross validation of the reconstructions. The parameters of the CGNE scheme 
are as before. Moreover, we show the residuals of the CGNR scheme for the 
least squares problem (1.2) with A^ = 256 (underdetermined) and N = 64 
(overdetermined) . 



19 



6 Conclusion 



Wc have shown that the optimal trigonometric interpolation problem at g-separated 
nodes in d dimensions is well conditioned for a polynomial degree N > 2dq^^ . However, 
in our further extensive numerical examples we observe that for A'^ ~ M^/"^ one can expect 
fast convergence of Algorithm 1. If we assume furthermore a uniformity condition q = 
cM~d for the sampling set X of cardinality M, then the total arithmetical complexity for 
solving the interpolation problem (1.3) up to a prescribed error is of order O (MlogM). 

We remark that dependent on the application, one solves the weighted approximation 
problem (1.2) or the optimal interpolation problem (1.3). Under some further mild 
conditions on the sampling set, both problems are solved efficiently by means of the 
conjugate gradient method in conjunction with the nonequispaced FFT. 
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